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Abstract 

The principles underlying plant development are extended to allow a more molecular mecha- 
nism to elaborate the schema by which ground cells differentiate into vascular cells. Biophysical 
considerations dictate that linear dynamics arc not sufficcnt to capture facilitated auxin trans- 
port (e.g., through PIN). We group those transport facilitators into a non-linear model under 
the assumption that they attempt to minimize certain differences of auxin concentration. This 
Constant Gradient Hypothesis greatly increases the descriptive power of our model to include 
complex dynamical behaviour. Specifically, we show how the early pattern of PINl expression 
appears in the embryo, how the leaf primordium emerges, how convergence points arise on the 
leaf margin, how the first loop is formed, and how the intricate pattern of PIN shifts during 
the early establishment of vein patterns in incipient leaves of Arabidopsis. Given our results, 
we submit that the model provides evidence that many of the salient structural characteristics 
that have been described at various stages of plant development can arise from the uniform 
application of a small number of abstract principles. 
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1 Introduction 

In the companion paper [14] we introduced a model for auxin dynamics that can predict key features 
of leaf vein patterns. While it is remarkable that diffusion of the hormone provides a sufficient basis 
to describe these patterns, this transport mechanism does not accurately describe the flow of auxin 
in plants. Facilitated transport of the hormone has been postulated since the proposal of the 
chemiosmotic theory, and recent molecular genetic and cell biological findings support this theory. 
Genes and related proteins have recently been identified as either influx (e.g., AUX) or efflux (e.g., 
PIN) auxin 'carriers'. Most interestingly, the experiments suggest that the hormone may in fact 
be 'pushed' against a concentration gradient — that is, opposite to the direction of diffusion flow 
in the sense of our earlier model. The carriers also appear to increase flow in the direction of 
diffusion flow at times, but why this dual function exists and when it shifts from one to the other is 
not well understood. Opinion among experimentalists is split, resulting in a conundrum for plant 
developmental biology. Recent experiments that focus on facilitated transport due to the PINl 
gene provide another dimension to the conundrum. Scarpella et al. [51] visualize PINl expression 
domain, or PED, at different stages of early leaf development and observe that eventually the PED 
defines the vein pattern exactly. Curiously, the PED undergoes several transformations whereby 
certain parts appear temporarily only to disappear just before the final pattern is established. Yet, 
the patterns that the PED defines are robust to cell divisions. The dynamics are non-linear and 
sometimes counter-intuitive. It would appear, therefore, that a carefully timed complex genetic 
machinery must exist for the purpose of controlling these events. 

Contrary to this intuition, in this paper we demonstrate that the underlying principles that 
guide such intricate events need not be many. While our focus remains abstract, our goal is to 
bring our model closer to cell and molecular biology. We concentrate on the phenomenon of polar 
transport rather than its molecular implementation and extend our earlier model to accommodate 
this abstraction. Building on the two principles from our earlier work, the Constant Production 
Hypothesis and the Proportional Destruction Hypothesis, we develop the Constant Gradient Hy- 
pothesis, which we use to derive a new computational tool to test the theoretical predictions. We 
introduce a non-linear model mathematically and then use it to analyze the observations reported 
by Scarpella et al. [51]. We focus on the shifting patterns of the PED during early leaf devel- 
opment, and discuss schematic simulations to offer an explanation for key stages of those events. 
Thus we illustrate how our theory captures certain intricate patterning events. Insofar as the polar 
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Figure 1: Sketch of predictions. 
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Table 1: The variety of auxin facilitators reveals enormous potential for complex interactions. Since 
many of these remain to be studied, we study their net effect. 

transport mechanism is implemented by a single substance, our simulations predict the behavior of 
that substance; however, (and this is key) our polar transport extension is meant as a simplifying 
abstraction of the Chemiosmotic Theory that captures the total effect of various substances acting 
together, rather than the behavior of any one substance. To preview the scope of our predictions, we 
sketch in Fig. 1 certain of the pivotal events in vascular development. Once our theory is developed, 
we return to these predictions via simulation. 

2 Auxin Transport 

Auxin travels through plants by moving between cells. It may leave a cell interior, the cytoplasm, 
pass through the plasmalemma and then reach the cell walls, or apoplast. It can then move into the 
cytoplasm of an adjacent cell and continue diffusing in this fashion throughout the plant. Perhaps 
the first real breakthrough in understanding how this movement occurs was the introduction of the 
chemiosmotic theory in the 1970s. Based only on physical chemistry and a handful of measurements, 
Rubery and Sheldrake [46, 47] and Raven [44] were able to predict that auxin cannot travel by 
diffusion alone; that some sort of active transport (i.e. one that requires energy expenditure) must 
be present; and that some sort of carrier substances must exist that move auxin in a preferred 
direction — toward the cytoplasm (influx) or away from it (efflux). 

It is "remarkable how accurately [this] molecular model [...] fits with the recent molecular 
genetic and cell biological findings" [58]. Increasingly, the evidence suggests that auxin transport 
is facilitated in both directions. The AUX proteins are a family of putative influx carriers [2, 56, 
13, 60, 30], while the PIN family [30, 58] and the multidrug resistance/p-glyvoprotein (MDR/PGP) 
family of proteins [7] are thought to be efflux carriers. Recent molecular techniques, which can 
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effectively 'color' specific molecules with sub-cellular precision [19], have shown that both types of 
facilitators may localize asymmetrically on the plasmalemma although MDR/PGP proteins do not 
always do that [50]. Indeed, many proteins have been implicated in auxin transport, as Table 1 
summarizes, and often more than one substance is present in the same cell at the same time. It 
is therefore difficult to infer the pattern of auxin transport from this detailed mechanistic view of 
auxin transport. 

In order to analyze how the patterns of auxin transport form we shall consider polar transport 
in a cell as the total effect of polar transport as predicted by the chemiosmotic theory. So, for 
example, if both AUX proteins and PIN proteins are expressed in a cell in equal measure and on 
the same parts of the membrane, then the action of the influx carriers will be canceled by the action 
of the efflux carriers and the total effect will be nil. Conversely, if we observe a total efflux effect 
in a given cell, then this only implies that the balance between influx and efflux carriers favors 
the latter. This simplification requires a refinement to our earlier model, which we prove to be 
consistent with the Chemiosmotic theory in Appendix A. In particular, we may limit the analysis 
to the movement of auxin between cell interiors and yet capture a sufflciently accurate picture of 
the whole process. We note that, if the molecular carriers in a given region of cell are always of 
the same kind, as seems to be the case in Scarpella et al. [51] (all PINl), then any predictions 
using this simplification apply directly to that carrier. For this reason and because most of the 
experiments that we shall analyze involve only PIN proteins, we shall refer to our simplification of 
directed (polar) transport as "PIN" transport. 

The key difference between the model that we develop here and our earlier formulation is that 
polar transport may move auxin against concentration gradients (from low to high concentrations) 
and therefore counteract the passive diffusion flow. The modeling challenge thus moves beyond that 
of Fickian diffusion, but its importance increases dramatically. The reward is that the capability 
for pattern formation is enlarged as well. We now explore this increased capability. 



We now consider some of the experimental evidence concerning PINl. This will be the basis for 
our assumptions about polar transport recognizing the fact that the experiments reviewed below 
show PINl but other carriers may be present in the same cells. 
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3.1 A Conundrum 

The appearance of PINl in detectable quantities is intimately related to auxin [59, 41], but the 
polarity of the carrier need not follow the diffusive auxin flow. For example, Reinhardt et al. [45] 
report that the pattern of PIN expression changes after a micro-application of atixin in shoot 
meristems — the polar transport acts toward the application site. In leaves, on the other hand, 
Scarpella et al. [52] applied the hormone to a small portion of a developing primordium and noticed 
that PINl oriented so that the hormone would be better drained away from the leaf, i.e. away from 
the application site. The application site in both cases has a higher auxin concentration than the 
rest of the organ so, considering each experiment separately, two mechanisms come to mind: 

Mechanism 1 : New PINl appears at a membrane against diffusion flow when this flow, or differ- 
ence of concentration, is large enough. 

Mechanism 2: New PINl appears at a membrane 'helping' diffusion flow. 

At first glance it looks like neither rule can explain both experiments because the phenomena 
appear to be mutually exclusive. But it turns out that one of these schematic mechanisms is the 

foundation of a model that predicts both phenomena in detail as well as vein patterns and the 
growth of organs more generally. 



Our first clue comes from considering what happens when a cell with PINl expression divides. 
In leaf primordia, the cells that will eventually become part of the midvein are marked by PINl 
expression that localizes toward the base of the cell, i.e. toward the stem and away from the leaf 
tip. This pattern appears as soon as the primordium emerges from the shoot apical nieristem and 
the polarity of the cells is maintained even during frequent cell divisions. When a cell with PINl 
expression divides perpendicularly to the midvein strand, a new membrane is created that also 
acquires PINl expression with the polarity of the original cell. Mechanism 1 fails to predict such 
behavior, while Mechanism 2 maintains the proper PINl polarity. We reason as follows. 

Consider the illustration in Fig. 2A-D. Initially, the concentration of auxin in the middle cells 
is higher than in the cell on the left and lower than in the cell on the right. When the middle cell 
divides its membrane is 'split' and a new barrier forms. The existing auxin carriers then drain auxin 
away from the daughter compartment on the right and push auxin into the neighbor cell on the 
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Figure 2: Illustration of the effect of two presumptive mechanisms for PINl creation in the context 
of a cell division. Three cells are shown in A and E with consistent polarities indicated by open 
arrows. When the middle cell divides, its membrane remains the same except for the new membrane 
that splits the cell in two daughter compartments. We take this to imply that the existing PIN 
structure is maintained and the problem is to determine how PIN is established along the new 
section of membrane. The concentration of auxin (denoted [lAA]) is plotted above the depiction 
of each cell. PINl, shown in green, is created according to one of two mechanisms: (1) toward 
higher concentrations (see C and G), or (2) toward lower concentrations (see D and H). Note that 
Mechanism 1 breaks the continuity of PIN polarity in both cases, whereas Mechanism 2 maintains 
this continuity. (A— D) PINl acts against diffusive flow, toward increasing concentration, e.g. as 
in the Arabidopsis root tip. (E— H) PINl 'helps' diffusive flow, e.g. as in leaf primordia. (A,E) 
Configuration before cell division. (B,F) Configuration during cell division (dashed line) but before 
PINl is established at the new interface. (C,D,G,H) Configuration in middle cell (and respective 
daughter cells) after PINl is created according to Mechanism 1 (C, G) or Mechanism 2 {D,H). 

left. On the other hand, the carriers from within the neighbor cell on the left push auxin into the 
daughter compartment on the left. Thus, the left compartment acquires a higher concentration than 
the right compartment. If PINl were to be created toward a higher concentration — Mechanism 1 
and Fig. 2C — then it should form in the right compartment thereby creating a bipolar cell. This 
would be the outcome after most cell divisions and, in particular, in the midvein strand where no 
bipolar cells have been found. By contrast, if PINl were to be created predominantly toward lower 
concentrations of auxin — Mechanism 2 and Fig. 2D — then the continuity would be maintained. We 
therefore conclude that the first appearance of PINl in cells is most likely due to Mechanism 2. 

Recall Schema 1 from [14] which postulates that a ground cell becomes c-vascular whenever it 
measures a sufficiently large Ac. At the schema level of abstraction, the change is realized as an 
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increased diffusion coefficient through the relevant interface, and the Ac there decreases. This is 
precisely the effect of Mechanism 2! In other words, we have now shown that Schema 1 can be 
implemented by auxin carriers. Therefore our derivation of Schema 1 informs the appearance of 
auxin carriers as a result of organ geometry and cell size distribution. The agreement between 
these predictions and experimental fact reinforces confidence that Schema 1 is implemented by 
these carriers. 



4 The Constant Gradient Hypothesis 

Even though Mechanism 2 will be a useful guideline to understand PINl dynamics, it does not 
explain how the carrier concentration at the cell membrane is maintained. Indeed, the protein cycles 
between the interior of the cell and the plasmalemma [41], and it is somehow destroyed during the 
early stages of leaf development [52] and when it undergoes reorientation [45]. Moreover, even if 
the first appearance of PINl is in the direction of decreasing auxin concentration, there is carrier 
presence against the auxin gradient. For example, the polarity of provascular cells near the distal 
tip of the Arabidopsis root is toward the tip, and the auxin concentration increases in the same 
direction [30]. This suggests that the mechanism responsible for PINl density at the plasmalemma 
exhibits some sort of hysteresis memory. Taken together, the cycling, the plasticity of expression, 
and the memory properties of PINl imply that the mechanism is a dynamical system governed by 
a differential equation. 

Mathematical formulations of this sort are common and find many applications in biology [43] . 
For example, they are the theoretical basis for genetic switches [3] where the equations can be 
understood in terms of two competing processes: the activity of promoters described by a production 
function, and the activity of inhibitors and natural degradation described by a destruction function. 
The production function depends on the concentration of the gene product whenever the product 
promotes the expression of the gene, and it can usually be described by a Hill function. If there are 
several promoter sites and cooperation between them, then the Hill coefficient is large [25, 23, 33, 17] . 
On the other hand, the destruction is a linear function of concentration whenever no inhibitors 
exist and only natural degradation depletes the substance. Configurations of this kind exhibit the 
behavior of a switch: if an external event lowers the concentration below a threshold r, then the 
system maintains the concentration at level A; if another external event pushes the concentration 
above this threshold then the system maintains the concentration at level B {B > t > A). If A is 
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biologically negligible (insufficient for additional reactions) but B is not, then we can say that the 
gene is turned 'off' or 'on' depending on whether the product concentration is below or above the 
threshold, respectively. 

There is more direct and independent support for the dynamics of PINl to be well described by 
a Hill function and by the concept of cooperation. Cooperation is central in the empirically derived 
Canalization Hypothesis by Sachs [49] whereby "flux begets flux." If the increase in flux is the result 
of "PIN" , then the carriers must exhibit an autocatalytic behavior because they are saturable: they 
have a maximum capacity for transporting auxin through the plasmalemma [42, 58]. Note, however, 
that we are referring to the net effect and that the actual behavior may involve more than a single 
substance. For example, it may be the result of mutual cooperation between the PIN proteins and 
the MDR/PGP proteins [37, 42]. Therefore to increase the flux of auxin due to those facilitators, 
the concentration of the molecule must be increased. We submit the following hypothesis: 

Hypothesis 1 (PINl Cooperation). PIN carriers at the plasmalemma attract additional PIN car- 
riers near them. The protein leaves the membrane when it has transported its load of auxin. 

In other words, the density of PIN carriers at a patch of the cell plasma membrane is determined 
by a dynamical system with a production function that depends on that density, a Hill function. 
The destruction function also depends on the density but can be assumed to be linear, because only 
a portion of the loaded proteins — those immediately adjacent to the membrane, those attached to 
it — can facilitate hormone transport at any given time. In effect, this configuration constitutes an 
auxin carrier 'switch.' 

As we argued above, the external events that turn the switch 'on' or 'off' depend on auxin 
and, particularly, on the difference of auxin concentration Ac through an interface. The rule which 
we called Mechanism 2 above implies that a large Ac in the direction of diffusion — i.e. a positive 
Ac — should increase the production function. On the other hand, the reorientation of the carrier 
following external application of auxin implies that a high gradient against the action of PIN should 
increase the destruction function. Yet, the strong expression of PINl against the auxin gradient in 
provascular cells of the root tip implies that the destruction function overwhelms the production 
function only when Ac is sufficiently negative. We therefore submit that PIN dynamics are modified 
by auxin in the following fashion: "'^ 



-"^In fact, a more mechanistic hypothesis that justifies Hypothesis 2 can be made. It is about both influx (e.g., 
AUX) and efflux (e.g., PIN) carriers and requires them to function in a dual fashion: Ac makes it difficult for PIN 
to unload auxin and easy for AUX. Here, Ac is the more immediately available difference of auxin concentration 



Draft of May 28, 2009 [14:28] 



5 COMPUTATIONAL MODEL OF ACTIVE TRANSPORT 

Hypothesis 2 (Effect of Aiixin on PINl Dynamics). PIN accumulation on the plasmalemma is 
promoted by Ac regardless of sign. The destruction function increases only when the diffusive flow 
opposite to the action of PINl increases, i.e. if the carrier pushes auxin against diffusion. 

In other words, PINl is more likely to be destroyed when it pushes auxin against a concentration 
gradient. But even in that case, the production process may still be sufficiently strong to compensate 
for the action of the destruction process and maintain carrier presence or, as in the root tip, even 
increase it. However, if the gradient is too steep, then destruction overwhelms production and 
carrier presence disappears. 

The lack of experimental data prevents us from determining exactly when this happens and 
the explicit forms of the functions that describe the dynamics. This underlines the advantage of 
our abstraction, though: we can use our theory to infer a functional role for auxin carriers that is 
consistent with our hypotheses and that will allow us to predict the distribution of polar transport 
intensity in plants. Some data of this sort are available, and more precise measurements can be 
obtained with current techniques. 

Recall that we required veins to possess an efficient transport mechanism for auxin. In particular, 
no concentration peaks should form there. On the other hand, we saw that the appearance of auxin 
carriers at a cell interface is most likely triggered by a high Ac and, initially at least, helps diffusion. 
In effect, carriers decrease this Ac. Working against diffusion, too, they maintain a Ac that does 
not exceed some fixed threshold because otherwise carrier presence would disappear. Therefore, we 
propose a third uniformity principle: 

Hypothesis 3 (Constant Gradient). Facilitated transport maintains the difference in auxin con- 
centration, Ac, constant. 

5 Computational Model of Active Transport 

The Constant Gradient Hypothesis suggests an elaboration of the computational model developed 
in the previous paper [14]. In this section we derive the mathematical and computational tools that 
are necessary to evaluate the theory. We develop a framework for simulations, which we use in the 

next section to illustrate how to explain experimental evidence of certain patterning phenomena. 

between the cytoplasm and the cell wall instead of, as we conjecture in Hypothesis 2, a Ac between the cytoplasms 
of two neighboring cells. 
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The tools that we develop here have a larger scope than is shown in the next section, as can be 
seen in the next paper [15] of this series. 

5.1 Background 

Even though axixin transport may be due to energy expenditure, the net effect may still be well 
described by Pick's law. This is the case when carriers are located in a non-polar fashion (homoge- 
neously) along the cell membrane. The PGP family appears to exhibit this behavior so it can be 
seen as a mechanism for improving the effective diffusion coefficient. In this case, our Helmholtz 
formulation from [14] applies and large domains of cells can be described by the continuous equation 

-=DV^c+--ac. (1) 

Thinking of cells as individual compartments, this means that the concentration c{i) of auxin in 
each cell i changes as a function of time according to the transport to neighboring cells Nbr(i), the 
production K/S{i), and the destruction —ac{i); in symbols: 

^= Dic{j)-c{i))+^-ac{i). (2) 



Diffusion Transport 



On the other hand, the PIN family of proteins (and PINl in particular) localize asymmetrically 
along the membrane — with negligible or no expression on one side of the cell and strong expression 
on the other. This is the scenario that we now study because it yields dynamics for auxin transport 
that deviates from Pick's Law: circumstances exist in which it can 'push' the hormone against 
concentration gradients. A typical assumption regarding the effect of transport facilitators such 
as PIN is that they add a linear term to the model. This turns out to be consistent with the 
chemiosmotic theory, as we show in Appendix A, so we incorporate this idea into our Helmholtz 
Model of [14] as follows: 

^-^= D{c{j)-c{z))+ J2 pT<j)- J2 pT<i)+^^-^<^- (3) 

j&Nbr{i) jeNbr jeNbr ^ ' 

^ ' ^ V ' 

Diffusion Transport Polar Transport 

Here p'" is the (active) polar transport coefficient for auxin being 'pushed' into the current cell i 
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from a neighboring cell j, while is the polar transport coefficient for auxin being 'exported' 
from the current cell toward the neighbor j. Since these coefficients represent the asymmetrically 
located putative carrier PINl, they also represent a measure of the work that the carrier performs. 
In particular, if PINl works at saturation, then the p's are proportional to PINl concentration at 

an interface. Thus, the biological appearance of PINl translates into the appearance of non-zero 
p's in this model. But how can we predict when that happens? The Constant Gradient Hypothesis 
offers an answer, provided that we can solve a non-linear transport problem. 

5.2 Formulation of Non-linectr Transport Problem 

Suppose the cell sizes in the domain of interest (e.g., a young leaf) are known and that we also 
know which cells express PINl. We are given the interfaces where the carrier is present and the 
direction in which it facilitates auxin transport, but we do not know how much protein is present 
there. ^ PINl may work either in the direction of diffusive flow or against it: these are two distinct 
modes of operation. We shall see how to assign an operation mode to each interface in the following 
section. Here we discuss a mathematical tool that allows us to compute the auxin concentration at 
equilibrium, c, assuming the modes of operation are known everywhere. Finally, this c is used in 
conjunction with Eq. 3 to estimate the amount of PINl present at each interface. 

In the first mode of operation PINl maintains Ac < n because it works with diffusion. We 
define the following non-linear transport function; 



!!>"'(Ac) 



where Dfast ^ D. Note that the usual transport function due to Fickian diffusion is (/)^*'^'^(Ac) = 
DAc. The only difference is that the diffusion coefficient increases non-linearly when the difference 
in auxin concentration through the interface exceeds the threshold ti. The dynamics of cell i 
become: 

E -^Wi) - cW) + - ac(i) 

ji£Nbr{i) ^ ' 

where Nbr{i) denotes the neighboring cells of cell i, and the fiux <j) now depends on whether the 
interface contains PINl helping diffusion or not, i.e. it is either 0"' or The effect of this 




^Note that our previous model [14] can make this sort of prediction. Also, the distribution of PINl can be obtained 
from experiments such as those in [52]. 
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formulation is that an interface with transport function (describing the flux of auxin) 0"' wiU have 
a Ac < Ti or Ac w n at equiUbrium for a sufliciently large -D/ast- This is equivalent to what the 
Constant Gradient Hypothesis postulates for PINl working with diffusion.^ 

When the carrier works against the concentration gradient, however, the hypothesis implies that 

the |Ac| is maintained below but near T2. For the purposes of the present paper, this behavior can 
be achieved by a constant flux. For example, if PINl is present at the interface between cells i and 
j, facilitates transport of auxin from i to j, and c{i) < c{j) (i.e. PINl works against diffusion), 
then the flux due to PINl is a constant Puphui- In symbols. 



So, we only assume values for uphill polar transport and we can solve the non-linear system to 
obtain the distribution of auxin concentration at equilibrium (details in Section 5.3). Using this 
distribution we can obtain the polar transport contribution in all cells (details in Section 5.4) and 
obtain all the parameters of the polar transport dynamics at equilibrium (Eq. 3 with ^ = 0). In 
particular, we obtain the values of the polar transport coeflicients p'" and which are related 
to the density of PINl at the respective membranes. 

The above analysis, together with our Schema 1 fleshed out as the appearance of PINl, allows 
us to make much more detailed predictions. These predictions follow in Section 6 and Section 7. 

5.3 Solving the Non-linecir Dynamics Problem 

Now we show how to flnd the equilibrium concentration of auxin according to the non-linear trans- 
port formulation. We need to flnd c(i) for all cells i that satisfles: 



where Nbr{i) denotes the neighboring cells of cell i; (j) is the flux of auxin from cell j into cell i; 
p{i) = is the production function; and a is the destruction constant. There are three types 
of flux functions of interest: diffusion according to Pick's Law c^^*^*^; non-linear diffusion 0"'; and 
polar transport ^"''"'^ {c{j) , c{i)) = ±PuphiH where the sign is positive if PINl facilitates transport 

^Note, however, that the non-linear dynamics do not say anything about how the system behaves outside of 
equiUbrium. This is simply meant as tool to compute the equilibrium. 





(4) 



jZNbr{i) 
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uphill from cell j into cell i and negative in the opposite direction. All three share properties that 
guarantee the dynamical system to have a unique solution. 

Proposition 1. Let (j){c{i),c{j)) have the following properties: 

1. Antisymmetry: 4'{x,y) = —<p{y,x) 

2. Monotonicity: x — y <v — w (j){x, y) < <f){v, w) 

3. </'.= f^>0. 

Then, there is a unique c that satisfies Eq. 4 which can be found using a Newton method. 
Proof. Let c be the distribution of concentration and define the dynamics as: 



jeNbr{i) 



Now the Jacobian with respect to c{i) is a matrix 

'd{ctm 



j = 



dcij) 

with off-diagonal entries, i ^ j, Jij = ^'^^'gc{j)^^^^ ~ ?^x(c(i),c(i)). The diagonal entries are 

T d{ct){i:) ^ d<f){c{j) , c{i)) \^ A (^(^\ ^ 

^ ' jeNbr(z) ^ ' jeNbr{i) 

where (f)y denotes the partial derivative with respect to the second argument. 

Now, using the antisymmetry property, differentiating with respect to x on both sides, we see 
that (j)x = —(f>y Therefore, the diagonal entries of the Jacobian become 

Jii = -0x(c(j),c(i)) - a . 

jeNbr{i) 

Thus, the third property implies that the diagonal entries are strictly negative and that they 
dominate the off-diagonal entries (because J^jjU'^ij — SjeAr6r(i) ?^a:(c(i), c(i))). As a result, the 
Jacobian J is always invertible and Newton's method applies. In particular, this means that 
J(ct) = (ct) = so such a c exists because the c > and ^^c < ^^pjoL under these 

dynamics. 
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It is also unique. Let the transport operator $ be defined as 

jeNbr{i) 

and suppose that both x and y are solutions. Thus, 

— ax = —p = — ay ^x — ^y = a{x — y) 

So the inner product 

{<^x — ^y, x — y) = a{x — y,x — y) > . (5) 

On the other hand 

{^x -^y,x-y) = {^x, x) - {^x, y) - {^y, x) + y) (6) 

where 

{<^x,x) = 'Yxi{<^x){i) ='Yxi i ^ (j){xj,Xi) 

i i \jENbr{i) 

Each pair of neighbors appears exactly twice in the flux function: once for each direction. Thus, 
the term i) corresponds to the term Xj4>{xi,Xj). Therefore, using the first property we can 

pair the two terms as {xi — Xj)(j){xj,Xi). The expression for the {^x,x) then becomes 

{^X, X) = ^ {Xi - Xj)(j){Xj, Xi) 
j^i Nbrs 

where the sum is over all pairs of neighbors i and j. Using the same argument, we obtain expressions 
for the other inner products and rewrite Eq. 6 as 

{^x — ^y,x — y) 

= Y.{xi - Xj)(i){xj,Xi) - {yi - yj)(i){xj,Xi) - {xi - Xj)(t){yj,yi) + {yi - yj)(i){yj,yi) 
= E-(^i - Xj)(l){xi,Xj) + {yi - yj)(t){xi,Xj) + {xi - Xj)(l){yi,yj) - {yi - yj)(t){yi,yj) 

Let A = Xi — Xj, (t){A) = ^{xi, Xj), B = yi — yj, and <p{B) = <p{yi, yj). So we can rewrite the terms 
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in the summation as 

-A(t>{A) + BcPiA) + A(j){B) - B(I>{B) = -{A - B)(<^(A) - (j){B)) . 

The second property of <f) guarantees that {A — B) and <^(A) — (l){B) have the same sign, i.e. 
-{A - B){(I){A) - (l){B)) < 0. Therefore, {^x -^y,x-y)< 0. Thus, recalUng Eq. 5, we see that 
(i>a; — x — y) =0, which either means that x = y and we are done, or that ^x = ^y. In the 
latter case, we have 

$x — ax = $y — ay —ax = —ay x = y 
so the solution is unique. □ 

5.4 Computing the Polar Transport Contribution 

Now suppose that the distribution of auxin concentration at equilibrium, c, is known — e.g., com- 
puted as in Section 5.3 — and that the polar transport dynamics are given by 

dc K 

j^{i)= E {D{c{j)-c{i))--Ki^i + T,i^i) + — -ac{i) = Q (7) 

jeNbr{i) ^ ' 

where p{i) = K/S{i) and a are also known. Our task here is to compute the effective polar transport 
contribution: the iTi—,j. The unsaturated polar transport and the effective polar transport are 
related 

TTi^j = Pi^jC{l) (8) 

SO we can rewrite Eq. 7 in matrix form; 

Vc + Vc + p- ale = (9) 

where 2? is the diffusion matrix as described in [14]; V is the polar transport matrix; p is the 
production vector; and I is the identity matrix. The polar transport matrix is everywhere zero 
except for interfaces i\j with polar transport. Thus, if the direction of transport is i — > j, then 
'Pij = Pi^j and Vii = —Pi^j + Other (where Other refers to other polar interfaces that cell i may 
have). This contributes —pi^jc{i) to the dynamics of cell i and pi^jc{i) to the dynamics of cell j as 
required. The unknowns are the Pi_>j, i.e. the entries of V. Let p be the vector of those unknowns 
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and write 

{V*)p = Vc . (10) 

Eq. 9 implies that 

pc = {-V + al)c - p, (11) 

which contains no unknowns and can be computed. The matrix V* is the dual of V and can be 
written as follows. The i*'^ coordinate of the vector Vc consists of the sum of the contribution of 
all outgoing polar transport activity — cell i exporting to neighbors — and all incoming activity — 
neighbors exporting toward cell i. 

jEout jEin 

SO combining with Eq. 8 we obtain an expression iTi^j 

jEout jEin 

where 7? is the vector of effective polar transport variables, and 11 is the n x m incidence matrix 
for the polar transport graph where n is the number of cells and m is the number of edges. Thus, 
letting k denote the edge = (i ^ j), we have 

+ 1, if/=i 

Hifc = i -1, iil = i 

0, otherwise. 

The polar transport graph is a subgraph of the cell graph and consists only of those cells which 
are incident to an edge that corresponds to an unknown ni^j — there are n' < n such cells. In 
particular, only n' out of the n rows of 11 have non-zero entries. 

Proposition 2. We can compute tt — and therefore p — if and only if the polar transport graph is 
an undirected forest, i.e. m = n' — c (c is the number of connected components), and c was obtained 
from a transport rule cj) that satisfies j) = —(j>{j i) for each edge 

Proof. First, we note that the rank of 11 is n' — c by Proposition 4.3 in [9, p. 24]. Thus, the non-zero 

portion of 11 is square and invertible if and only \i m = n' — c. Suppose that 7? was computed using 
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this square matrix. The variables satisfy all n' equations, not only the n' — c. To see why, note that 
TTj^j appears in exactly two equations but with opposite signs. Also, diffusive transport works in 
the same way, as did the original transport rule. Thus, both sides of Eq. 11 satisfy the property 
and TT is a solution to the whole system Eq. 13. □ 

The exact nature of transport facilitation that PINl provides is unknown, but the effective polar 
transport provides a measure of the effect. Hence, assuming a model of auxin transport by PINl, 
we can compute how much protein is needed to achieve the effect through a relation analogous to 
Eq. 8. In the simplest model the density of PINl is only proportional to tt. 



6 Predictions of Vein Patterns in Leaves 

We shall now show that wc can not only predict where and when new c- vascular strands form, but 
also the relative strength of polar transport that implements the required improvement in auxin 
transport. The the best of our knowledge, PINl is a likely and, possibly, paradigmatic instance of 
molecular polar transport and a precursor to the venation pattern. However, the pattern of PINl 
expression, referred to as PED (or PINl expression domain), does not correspond exactly to the vein 
pattern. Recent experimental studies by Scarpella et al. [52] demonstrate that only a subpart of the 
PED becomes vein pattern. Surprisingly, some parts of the PED disappear. Here, we shall show 
that our theory is consistent with these observations by showing how each of the major empirically 
defined stages can be explained within our framework. We shall model the behavior of cells within 
a developing leaf — the so-called leaf primordium — and the adjoining cells of the shoot apex (see 
Fig. 3 and Fig. 4). We begin with a cartoon domain of cells, to show qualitative properties, and 
then proceed to compare these results to experimental data. 

6.1 Emergence of leaf primordium and formation of midvein 

Our model predicts the emergence of a midvein as a result of leaf growth. The formation of a new 
leaf primordium on the shoot apical meristem (SAM, see Fig. 3) begins by the reorientation of 
polarity in some cells and the formation of "convergence points" — red arrowheads in Fig. 5. This, 
as our model predicts, results in the creation of a locally extreme difference in auxin concentration. 
Fig. 5B illustrates that Ac is largest near the convergence point between the cell marked with a cross 
and its neighbor marked with a red star. The large Ac promotes PINl creation and, eventually, the 
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Figure 3: Structure of shoot apex and location of leaf primordium. (A) Arabidopsis thaliana 
seedling (from [53]). (B) A slice through the shoot apex of A. thaliana [31, Fig. 5] found under 
highlighted area in A. (C) Schematic of slice in B. SAM: shoot apical meristem, responsible for 
growth along the main axis. 




Figure 4: Setup for simulations. (A) Shoot apex of A. thaliana as in Fig. 3C. (B) Simulation 
domain corresponding to region of real organ. Red arrowhead points at tip of leaf primordium. 
Cell sizes increase moving away from the tip, compare with A where cells are largest in the stem. 
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B 













Figure 5: Predictions of auxin concentration during initial stages of primordium formation. (A) 
Magnification of the tip of our simulation domain. Carrier expression localizes asymmetrically and, 
along the marginal cells, points toward a convergence point — cell with red arrowhead. (B) The 
predicted distribution of auxin. Note that the largest Ac is between the cell marked with a red cross 
and its neighbor marked with a red star. This is the most likely place for new PINl. As A shows, 
these cells are at the tip of the inward forming midvein PINl expression domain. Exaggerated c 
for display purposes. 




I 



ts 



Figure 6: Prediction of midvein from MPED using Algorithm 1. The domain and the marginal 
PINl expression, green bars, are shown. Also shown are four snapshots during the execution of 
Algorithm 1 labelled ^1,^2,^3, or ^4. The new PED forms where the midvein of the new leaf should 
be, green arrows. Note that this pattern is not prescribed, but the result of implementing our 
Schema 1 of [16] as Algorithm 1. 

interface on the side of the cell with a cross is endowed with the carrier. Following the extension 
of the middle PED, the extreme Ac moves to the next interface of the cell with the star and the 
procedure repeats. But each time the PED is extended in this fashion, the extreme Ac decreases 
until it is too small to create any new PINl (Section 3). 

To illustrate, consider Fig. 6. We are given the marginal PINl expression and follow the pre- 
dictions of our model through time as they give rise to the midvein. At time ti we show the 
configuration from Fig. 5. Then new PINl appears and creates a large Ac at the tip of the PED. 
Observe that this Ac is smaller at time ^2 than at time <i. It is even smaller at time and, at time 
t4, it becomes too small to create new PINl. So the midvein PED extends a finite distance from 
the primordium tip toward the stem below. 

In fact, these predictions are consistent with experimental observations. Our initial setup mimics 
the biology. Fig. 7, and our predictions for an extending midvein during leaf growth coincide with 
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Figure 7: Predictions of auxin concentration during initial stages of primordium formation. (A,B) 
Initial stages in the formation of leaf primordia (from [52]). Sub-cellular localization of PINl 
directed toward so-called "convergence points": cells shown with red arrowheads. (C) Model of 
this behavior. PINl, depicted in green, organize in marginal PEDs (MP ED) facilitating auxin flow 
toward the tip of the primordium (the convergence point) from all sides. The hormone is evacuated 
by the inwardly directed PINl expression that will eventually become part of the midvein in the 
new leaf. (D) Entire simulation domain. Selected region enlarged in C. (E) Distribution of auxin 
concentration at equilibrium computed as described in Section 5. The midvein PED works with 
diffusion — mode 1 — and the MPEDs work against diffusion — mode 2. The largest Ac is at the end 
of the midvein PED between the cell marked with a cross and the neighbor marked with a red star. 
This is the most likely place for new PINl. Exaggerated c for display purposes. 

the empirical patterns Fig. 8. But our predictions go even further. We compute the relative demand 
for polar transport along the whole domain under the Constant Gradient Hypothesis. Thus, if the 
concentration of PINl is proportional to this demand we predict the distribution of the carrier 
throughout the leaf. Fig. 9 provides one example. Note that more auxin carrier is needed in the 
midvein than in the marginal PEDs, which is consistent with PINl expression in a real leaf. 



6.2 Convergence points on leaf margin 

The pattern of PINl expression shown in Fig. 9 is maintained as the young leaf grows. In general, 
when a cell expressing PINl divides the daughter cells inherit the polarity of the mother cell. This 
rule is broken only occasionally, and when that happens a so-called convergence point results. Our 
model can explain this phenomenon in the following fashion. 

First, examine what happens when a cell on the margin divides. Fig. 10 illustrates the distribu- 
tion of auxin concentration after the new cell membrane (or wall) has formed sufficiently to present 
a barrier and the concentrations do not change. Notice, in Fig. lOB, that the largest difference 
in concentration is between the two daughter cells and that it suggests that new PINl should ap- 
pear so that the polarity is consistent with neighboring cells. This behavior is generic — Fig. 12 C 
demonstrates it for all cells along the margin but it also holds in the midvein. 

The situation in the midvein is somewhat different because there are more neighbors and the 
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Figure 8: Growth of leaf primordium induces growth of midvein. Note that as the midvein grows, 
it connects the tip of the leaf — where cells are small so c is high — to cells in the stem — which are 
large so c is low. Thus, the concentration at the tip of the leaf primordium decreases. In effect, the 
midvein acts as a straight bar being pushed up at the tip but pulled down near the stem. This is due 
to PINl following Constant Gradient Hypothesis. Three primordium sizes are shown in increasing 
order: A~D, E~G, and H-K. (A,B,E,H,I) from [52]. 




Figure 9: Prediction of PINl strength computed as described in Section 5. (A) Sample leaf from 
[52]. PINl expression in green, stronger just under tip. (B) PINl polarity shown as green bars. 
Flow due to polar transport shown as green arrows. (C) Prediction of PINl concentration. Height 
of green bars denotes concentration. The expression is low along the margin and high in the midvein 
for both the real leaf, in A, and predictions, in C. 
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Figure 10: Auxin distribution following a cell division on the leaf margin. (A) The domain. 
Rectangular area blown up in B. (B) Ac shown as a yellow segment starting from the cell with 
the higher concentration. E.g., the largest Ac = c(2) — c(l). Since PINl is created due to large 
differences in concentration and in the direction of that concentration, the new PINl accumulates 
inside cell 2 toward cell 1. Polarity is maintained consistent with neighboring cells. (C) Auxin 
concentration shown as a height map over the domain. 



strength of PINl is largest in the midvein. Still, Fig. 11 shows that PINl polarity can be maintained 
in the midvein for much the same reasons as before. However, the ensuing difference in concentration 
through the newly formed interface is significantly larger for a dividing midvein cell than for a 
dividing margin cell. Fig. 12_D demonstrates this graphically. As a result, new PINl should appear 
quicker and have an effect on the concentration faster for dividing cells in the midvein than for 
dividing cells along the margin. Therefore, any PINl that may start forming at interfaces other 
than the new one would probably not form in sufficient quantity to remain there after the new 
interface has acquired its share of PINl.^ The cells along the leaf margin, on the other hand, 
cannot create Ac quite as large and the ensuing slower acquisition of PINl can result in PINl 
appearing — and persisting — at interfaces other than the newly created one. 

The most interesting phenomenon of this sort is the flip of polarity at a nearby interface, which 
creates a convergence point: a stable discontinuity of PINl polarity along the margin. Fig. 125 
illustrates the differences in concentration resulting from a cell division that created cells 2 and 3. 
Besides the largest Ac through the new interface, there are other interfaces with non-trivial Ac. 
Fig. 12 C plots three of the largest Ac as a function of the location of the dividing cell along the 
margin. Of the three, only the Ac through interface 1 — > 2 can break the continuity of polar cells 
and create a "flip" . The most likely location for such a flip is where Ac is largest and it is possible 



^In other words, the fast action of PIN at new interface will not allow the Ac at neighboring interfaces to exceed 
ri long enough to make p > p-r at the neighboring interfaces. Here ri is such that the PINl dynamics can, given 



sufficient time, switch to the "on" position, i.e. [PINl] > pr- 
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Figure 11: Auxin distribution following a cell division at the midvein PED. (A) The domain. 
Rectangular area blown up in B. (B,C) Largest Ac through new interface as in Fig. 10. 



for a flip to remain there. Fig. 13 illustrates a likely sequence of events that results in a stable flip of 
polarity. In particular, this gives rise to the so-called convergence point on the margin. The stability 
of this discontinuity is due the relatively slower dynamics of PINl along the margin as compared 
to the midvein are important here. If a flip were to occur in the midvein, then the resulting Ac 
against the action of PINl will grow faster than the neighboring interfaces can compensate and, 
eventually, will make PINl flip back. Hence, the flip will be unstable in the midvein. 



6.3 Formation of the first loop 

Once a convergence point has been established the pattern of PIN may continue progressing inwardly 
in a manner similar to the midvein appearance (Section 6.1). On the other hand, the convergence 
point may remain stable for a while — i.e. observable in experiments — and, we predict, the strand 
connecting it to the midvein may emerge according to slightly different dynamics. Fig. lAA shows 
that the largest Ac is near the midvein while the peak of concentration is at the convergence point. 
Accordingly, the new strand will emerge at the midvein flrst and make its way to the convergence 
point — ^just like a new strand in an areole progresses from the existing vein toward the peak of 
concentration (e.g., see [16]). 

In fact, dynamics much like the ones discussed for areoles predict the next step in the elaboration 
of the flrst loop. The PIN expression domain in Fig. 14_B that includes the portion of the margin 
between the tip and the convergence point, the midvein and the newly created strand form, geomet- 
rically, an areole. Surprisingly, the distribution of auxin concentration in this geometric areole looks 
very much like the distributions from [16] and [14] as Fig. IAC,D illustrate: low concentration along 
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Figure 12: Cell division effects along leaf margin. (A) The rectangle encloses the cells used in 
C. Example of one cell dividing while all others do not. (B) The neighborhood of a cell that has 
recently divided. The daughter cells are labeled by the numbers 2 and 3, and the new interface 
between them does not yet express PIN. The curves in C are obtained by sliding this neighborhood 
along the margin (rectangle in A): think of the position of cell 1 as moving along the margin, while 
the other compartments have positions relative to that of cell 1. (C) Three curves obtained from 
measuring Ac as indicated. (D) Effect of cell division along the whole PED. Hormone distribution 
computed for each location separately assuming that a cell has divided at that location only (as in 
Fig. 10 or Fig. 11). Height of bar is Ac through the new interface without PIN. Aggregate plot. 
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Figure 13: Creation of a convergence point at the leaf margin due to a cell division. PINl is 
created where Ac is large in the direction of arrows; destroyed if against such arrows. (A) Initial 
configuration. Cell divided into cells 3 and 2. Expect new interface to acquire PINl 3 — > 2. 
However, Ac at 1 ^ 2 is large enough, so PINl will flip there. Also, 5^2 will be created. (B) 
The predicted configuration form A. The largest Ac is against the PINl acting 5 ^ 2 so polarity 
will flip. (C) Updated configuration from B. (C,D) If the supporting PINl can appear fast enough, 
then the Ac through 2-^1 will not eliminate the original flip and the configuration will be stable 
as in E. (E) Stable configuration: there are no more high Ac. Cell 1 is a bipolar cell because PIN 
expresses on both its left and right interfaces. Cell 2 is a convergence point because cells on both 
sides exhibit polar transport toward cell 2. 

boundary, peak inside, and highest Ac at the boundary. The new strands inside this geometric 
areole can be predicted by following the largest Ac as discussed in [16] and the companion paper 
[14] — we obtain Fig. 15. 

The key difference, however, is that the boundary PED of the areole is not entirely stable. 
The interface of the bipolar cell with the convergence point actually expresses PIN in the direction 
against diffusion. Thus, a sufficiently quick infusion of auxin into the convergence point cell can 
increase Ac enough to cause a fiip, remove the bipolar cell, and reestablish the continuity of PINl 
polarity along the margin. Such an infusion can result from the creation of the new strand. 

Indeed, as Fig. 15 illustrates, the extension of the PED in the geometric areole can remove 
PIN from some interfaces. As the new PIN expression emerges (Fig. 15^), the right hand side 
strand is infused with additional auxin. Even though the strand may potentially adapt to the 
new demands — by increasing the amounts of PIN that drain the hormone away — any adaptation 
requires time during which the whole strand will see higher overall concentrations. This is especially 
important at the convergence point where, in particular, the Ac through the interface with the 
bipolar cell (red bar in Fig. 15) will rise. Since, as we argued in [14], the new strand (solid arrow) is 
likely to extend in quick bursts, the infusion of auxin into the CP is likely to increase this Ac and 
fiip the PIN expression. Once that happens, the PED connecting the margin to the midvein will 
gradually disappear as well — the stability of the convergence point requires all of the supporting 
structure as illustrated in Fig. 13. 
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Figure 14: Emergence of strand connecting convergence point (CP) to niidvein; creation of geometric 
areole. (A) Distribution of Ac. Notice that the Ac at the convergence point is zero so there is a 
peak of concentration there. Conversely, the largest Ac is at the midvein, so a new strand could 
emerge from the midvein and end at the convergence point by following the usual local rules. BC: 
bipolar cell; CP: convergence point. (B) Domain with new strand. Rectangle encloses the geometric 
areole. (C) Close up of geometric areole. Note the red bar under the interface between the BC and 
the CP. (D) Concentration, c, of auxin in geometric areole. Red bar corresponds to red bar in C. 
Note that c{CP) > c{BC) so PINl works against diffusion and an auxin infusion into CP can flip 
the polarity of the interface between BC and CP. 
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Figure 15: Formation of first loop. (A) The new PIN expression as predicted by our model is 
shown (solid arrow). This increases auxin supply into the strand along X, which causes the concen- 
tration to increase there. Since the midvein cannot drain the new infusion of auxin immediately — 
adaptation time — there is a back propagation effect (dashed arrow) which causes the concentration 
of auxin to increase even at the convergence point (CP). (B) The bipolar cell disappears since the 
PIN between BC and CP flips sides due to the increased Ac caused by the new strand. As this 
happens, the PIN at interface 1 starts pushing auxin against the gradient and eventually disappears. 
Next, interfaces 2 and 3 disappear for the same reason. (C) The new strand meets in the middle of 
the areole — where the concentration peak is — creating a bipolar cell there (BC) but the connection 
to the margin is now lost. 

7 PIN Patterning in the Arabidopsis Embryo 

The same Non-linear Polar Transport model which accounts for many key events during vein for- 
mation also explains key patterning events at the earliest stages of plant development: the embryo 
formation. Fig. 16 shows how to obtain the geometry of embryos from images of the tissues. Then, 
in Fig. 17 and Fig. 18 we show how the polar pattern of PIN expression is predicted by our model. 

We begin with the assumption that no polar pattern exists. Uniform expression of PIN is 
equivalent to diffusion transport — as we argue in Appendix Section A — so our Hclmholtz model 
suffices to predict the initial distribution of auxin concentration for an embryo in the globular 
stage. Fig. 17 A, B. This distribution exhibits a maximal difference in concentration between the 
suspensor cells numbered 24 and 25 which, according to Schema 1, makes that interface the most 
likely first location for strong PIN expression. As this happens, the Ac there decreases and causes 
PIN to appear in several other suspensor interfaces, as shown in Fig. 17C. In turn, this event causes 
the Ac between the suspensor cell 24 and its neighboring embryo cells to increase significantly, and 
polarizes the expression of PIN at those interfaces. This polarization continues into the embryo — 
like the formation of a new strand in an areole — until it reaches the smallest cells: cell 5 and cell 
6 in Fig. 18^. At that point no significant Ac exists in the globular stage embryo. The triangular 
stage, resulting from several cell divisions, presents a different cell size distribution and significant 
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Figure 16: Obtaining the geometry of Arabidopsis embryos. (A) Globular stage (from [11]). (B) 
Domain extraction from A following the procedure in Section C. Embryo: cells 1 through 23; 
suspensor: cells 24 through 29. (C) Triangular stage (from [11]). 

Ac appear again, Fig. 18B. In fact, a convergence point appears near cell 11 due to the established 
pattern of PIN from the earlier stages. 

The final predicted pattern, shown in Fig. 18C and Fig. 19C, captures two important qualitative 
characteristics observed in the lab: one regarding the emergence of new leaves and one regarding 
the root tip. The theory developed here explains the initial formation of the 'convergence point' 
patterns that mark the appearance of new leaf primordia. In the case of the embryo, these are the 
cotyledons, but the same process could be taking place near the shoot meristem in mature plants. 
The theory also explains why the polarity of PIN in the root is towards the root tip. That is the 
result of the original differences of cell size. The suspensor cells, being larger, produce auxin at 
a lower per- volume rate than the embryo cells, so suspensor cells are predicted to have a lower 
concentration of auxin. Thus, the difference in concentration sets up the polar pattern near the 
root tip, which will be maintained through development due to the memory-like behavior of PIN 
dynamics. In a companion paper [15] we shall show that this initial pattern sets the stage for 
the observed 'reflux' pattern in mature roots, and that the Polar Transport model developed here 
explains the events that lead to the mature pattern as well. 
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Figure 17: Prediction of PIN pattern in embryos. (A) Concentration of auxin assuming Helmholtz 
conditions, i.e. no polar transport. (B) Differences of concentration, Ac, for auxin distribution 
as in A. Arrow size proportional to Ac. Note largest Ac between suspensor cells 24 and 25: PIN 
should appear there first. (C) Ac after introducing some of predicted polar transport (red arrows). 
Note large Ac between suspensor cell 24 and the embryo: PIN should form there next. 
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Figure 18: Prediction of PIN pattern formation in embryos (continued). (A) Predicted PIN pattern 
by following Ac. No more large Ac exist. (B) Ac in a triangular stage embryo with established 
PIN pattern. Note that the new cells (size distributions) create a non-trivial Ac between cells 20 
and 11: PIN should form there next. (C) Final PIN pattern predicted by our model. Note the 
'convergence point' at cell 11. Blue (color coded) arrows show Ac — size represents value. Red 
arrows represent polar transport — their size does not indicate strength. 
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Figure 19: Prediction of polar expression of PIN in Arabidopsis embryos and comparison to real 
measurements. (A) Globular stage. (B) Triangular stage. Note that the suspensor cells are larger 
than the embryo cells. Hence, our Helmholtz model predicts higher concentration in the embryo 
than in the suspensor. Black arrows point at where the root tip will be. Insets 1 from [11]; insets 
2 and 3 from [5]. (C) Our model predicts the typical polar transport patterns in both the leaf 
primordium and near the root tip: (1) the 'convergence point' pattern near the embryo tip where 
a leaf primordium will emerge (a cotyledon); and (2) the polarity of PIN expression is toward the 
root tip and due to the larger suspensor cells. 

8 Conclusion 

The diversity of phenomenology, data, and temporal information about plant development is dif- 
ficult to unify unless an abstract perspective is taken. In this series of two papers we sought to 
articulate several of the basic principles that could underlie the developmental biology of plants, and 
to place these within a mathematical framework that allowed their consequences and implications 
to be calculated. 

Realizing that distance information is fundamental, we introduced three principles that could 
work together within a reaction-diffusion model to predict how global distance information could 
be signalled locally, at the cell level. The first two principles, that auxin is produced at a constant 
rate in each cell and is destroyed in proportion to concentration, yielded models that predicted 
qualitative auxin distributions. Together with a simple schematic rule that cells begin to convert 
from a ground to a vascular state provided the concentration difference exceeds a limiting value led 
to a model of vascular formation. 

In this paper we elaborated this schema to a more realistic biophysical level. Facilitated trans- 
port of auxin by carriers such as those in the PIN and AUX families serves to effect the ground-to- 
vascular transition, but it does so in the context of non-linear dynamics with substantial predictive 
power. While it is widely accepted that PIN patterns prefigure vascular patterns, we were able to 
show that this dynamic — including the discrete and specialized structures that arise within it dur- 
ing plant development — can be predicted. Simulations showed how the leaf primordium emerges, 
how convergence points arise on the leaf margin, how the first loop is formed, how the intricate 
pattern of PIN shifts during the early establishment of vein patterns in incipient leaves of Ara- 
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Figure 20: Schematic depiction of cells (compare with [26, Fig. 1]). Five major compartments are 
shown: three cell interiors and two cell walls. The permeabilities through the plasmalemma are 
shown for all four compartment junctions and denoted by p{i). 

bidopsis, and how the midvein forms. Most importantly, we showed how the embrio provides an 
initial configuration sufficient to get it all started. 

Just as principles provide the foundation for theoretical physics, and mathematics provides 
the machinery to realize their consequences, we have sought to provide an abstract view of the 

development of pattern structure in plants. Our principles were designed to articulate key aspects 
of plant biology when viewed from an abstract, "high in the air" perspective. We believe they 
will hold as the facts relating to auxin production, destruction, and facilitated transport develop 
further; and that their predictions will lead to further experiments revealing the beauty and subtlety 
of plant development. 



A Simplification of the Chemiosmotic Model 

We now show that the net effect of active transport as predicted by the chemiosmotic theory may 
be captured by a small extension to our earlier model. Specifically, we show that the same effect 
may be obtained without the need to introduce separate compartments for cell interiors and cell 
walls — it suffices to keep the interior of cells as discrete compartments and to introduce the effect 
of carriers as directional active transport coefficients. These coefficients are a direct measure of the 
net effect of facilitated auxin transport, which we develop in the main body of the paper. 

Suppose, then, that we have a system of five compartments as in Fig. 20. Let J{j — > j + 1) 
denote the mass flux between compartments j and j + 1, so J{j — > j + 1) = —J{j + 1 — > j). Deflne 
$(j ^ j + 2) as 

Our task is to show that it is possible to think of $(j — > j + 2) only in terms of the compartments 
j and j + 2, the cell interiors. 
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We do this by examining the behavior of the system at equiUbrium. The substance is only 
produced inside cells and it is also depleted there by various metabolic processes. The following 
equations formalize these assumptions: 



(1) 


J{j- 


2 - 


- i - 1) + J{j - 
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(2) 




1 - 


^ j) + J{j + 1 - 
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(3) 


JU- 


* j ' 


f 1) + J0• + 2- 


-> J + 1) 
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(4) 


J{j- 


1 - 


■. J - 2) + p{j - 


2) - ac(j - 2) 


= 


(5) 


J{j + 


1 - 


i + 2) + p{j + 


2) - Qc(j + 2) 


= 



(15) 



The concentration of the substance in each compartment is assumed to be uniform and is denoted 
by, e.g., c{j) for the middle cell interior. Each equation corresponds to a compartment. The first 
three describe the cell wall c(j — 1), the middle cell interior c(j), and the second cell wall c(j + l). The 
last two equations describe the equilibrium dynamics of the leftmost and rightmost cell interiors, 
respectively. The production of the substance is captured by p and the destruction by ac with a 
a constant. Let T{j) = p{j) — ac{j) and write the dynamics (at equilibrium) of the cell interiors 
assuming that the cell membranes and cell walls together are a single interface and that $ describes 
the flux through those interfaces: 

(a) $(j-2^i) + $(j + 2^i)+ T{j) = 

(b) $(j^.?+2)+ T{j + 2) = (16) 

(c) cj>(j^j_2)+ T{j-2) = 

It is easy to check that Eq. 16 can be obtained from Eq. 15 by the following linear combinations: 

(a) .= (l) + (2) + (2)+(3) 

(c) := 

Thus, if the $ flows only depend on the concentrations in three of the flve compartments, then 
Eq. 16 has a unique solution which is compatible with the chemiosmotic theory. We now derive the 
conditions under which this holds and, therefore, demonstrate that our simulations are in agreement 
with the chemiosmotic theory. 

The substance is a weak acid and that the concentration c = [HA] + [A~], i.e. the total 
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concentration of protonated (associated) acid and the concentration in ionic form (dissociated). 
The dynamics of c depend on the dynamics of both [HA] and [A~], so we need to express each one 
as a function of c. Noting that the dissociation of a weak acid is given by the reaction 

HA^H+ + A-; K„=N^mol-r' 



we write 



which together with c = [HA] + [A~] gives 

1 IQpH-pK 

[HA] = ——rj — ^ — -c and [A~] = — — -c . (17) 

The change of total concentration due to substance transport is the result of two simultaneous 

effects: (1) the protonated form HA diffuses through the membrane, and (2) the anion A~ moves 
across the membrane as a result of auxin carriers and the electrical potential established by the 
difference in pH between the interior and exterior (cell wall) of the cell [26]. This dynamical 
behavior is therefore (see, e.g., [27]): 

J(i ^ e) = Pka {[BA], - [HA]^) + PA-g{V) {[A-], - [A-]if{V)) (18) 

where Pha is the permeability of the protonated acid; P4- is related to the amount of auxin carriers 
on the membrane; V is the membrane voltage; the subscript e = j + 1 and i = j denote the exterior 
and the interior of a cell, resp.; f{V) = exp(0) for ^ = —FV/RT, an electrical term consisting of 
the Faraday constant F, the gas constant R, and absolute temperature T; g{V) = i^fi^y) ■ We do 
not include the diffusion of A^ through the membrane because it appears to be negligible [4, 28, 29]. 

Let = 1/{10P"^-P^ + 1) and aj = 1 - Pj. Let h = g{V)f{V) and, using Eq. 17, rewrite 
Eq. 18 as 

J{3 ^ i + 1) = Pha {PAj) ' Pj+Aj + 1)) + PA{j)g Mf) ' faj+ic{j + 1)) 

It is reasonable to assume that fj = fj+2 because both functions relate the interior to the 
exterior of the cell in the same order and the voltage appears to be the same. However, the voltage 
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has opposite sign when the order is reversed, so that fj = —fj+i = l/fj+i. Therefore, setting 
f = fj, we obtain 

^ i+2) = Pha {Pj+ic{j + 1) - Pj+2c{j + 2))+PA{j+l)fg (aj+ic(i + 1) - jaj+2c{j + 2)^ 
We can now write 2$(j j + 2) = J{j ^ j + 1) + J{j + 1 ^ j + 2) as 

2$(i ^ i + 2) = Pha {Pjc{j) - Pj+2c{j + 2)) + P^(i)5a,-c(j) + {PaU + 1) - PaU)) fgaj+ic{j + 1) 

- PA{j + l)9aj+2c{j + 2) 

(19) 

Substituting PaU) = -Pa + we can rewrite Eq. 19 as 

2$(i ^ i + 2) = (/3jc(j) - f3j+2c{j + 2)) + PaQ (a,c(i) - a,+2c(j + 2)) 

+ PU)gajc{j) - p{j + l)gaj+2c{j + 2) + + 1) - fgaj+ic{j + 1) 

(20) 

Now if p(j + 1) = p(j) — in case both sides of the cell have the same permeabilities — ^then the 
last term becomes zero. In fact, even if the permeabilities are different the term can be ignored 

because fg ^ according to the best available estimates (i.e. V — — 120mV). Also, since the 
pH difference between the inside and the outside of a cell is assumed to be the same for all cells 
(metabolically maintained), we can write (}j = (3j+2 = /3 and a = 1 — (3. Eq. 20 becomes 

2$(i ^ i + 2) ~ Pha{1 - I3){c{j) - c{j + 2)) + PA9P{c{j) - c{j + 2)) + gl3{p{j)c{j) -p{j + l)c{j + 2)) 

so assuming w.l.o.g. that p{j + 1) = (i.e. because Pa = P{j + 1)) and collecting terms 

2$(i ^ i + 2) ~ (c(i) - c(i + 2)) {PHAii + PAgP) + gppUHj) (21) 

^ ^ ■' ^ — V — ■' 

Diffusion coefF. Polar transport 

Hence, the flux of auxin between cells can be written as the flux between cell interiors as 

2$(j ^ J + 2) ~ D{c{j) - c{j + 2)) + P{j)c{j) (22) 

Note that if the auxin carriers are homogeneously distributed on the cell membrane or if there are 
none, then P{i) = and the flux is governed by Pick's Law. If, however, the carriers are only 



Draft of May 28, 2009 [14:28] 



36 



C GEOMETRIC DOMAIN DEFINITION USING VORONOI DIAGRAMS 



located on one side of the cell, then = and P{i) in Eq. 22 describes the effect completely. 
In particular, P{i) is proportional to the amount of carriers on the portion of membrane of cell j 
adjacent to cell j + 2 (as in Fig. 20). 



B Midvein PED prediction 

Here is the algorithm used in Section 6.1. It is possible to predict the length of the midvein by 
measuring the length (or size) of the marginal PINl expression, the MPEDs. We fix a threshold 
Ti on Ac for the formation of new PINl and measure the MPEDs. We obtain c as in Section 5 by 
first assuming that PINl works against diffusion in MPEDs — mode 1 — and with diffusion in the 
midvein PED (if such a PED exists) — mode 2. We check each interface where PINl exists against 
the computed c and change its operation mode accordingly. If all interfaces are consistent with 
c, then we find the largest Ac in the domain and check whether it is above n. We extend the 
PED if Ac > n; otherwise, we finish the procedure. This procedure is collected in Algorithm 1 
and ilhistratcd in Fig. 6. Wc can predict the new PEDs in this way provided that the marginal 
expression of PINl is known. 



C Geometric Domain Definition using Voronoi Diagrams 

Assuming that auxin diffuses much faster inside the cell than through the cell wall, we do not 
need to model auxin transport inside the cell, only between neighboring cells. This is the basis 
of the simulations in Ref. [14], for example, but the geometry and topology of real cells is not 
accurately captured there. In particular, the neighbor relations are inaccurate as are the sizes of 
the interfaces. Our approach based on Voronoi diagrams solves this problem. We now describe how 
arbitrary geometric domains with the required properties can be defined from the Voronoi diagram 
of a collection of points, and then demonstrate how images of plant tissues can be segmented into 
cells using this approach. The pseudo-code for the procedure is presented in Algorithm 2.^ 

Consider Fig. 21A,B. Two dimensional points (open circles) are defined and the Voronoi diagram 
computed. The diagram consists of straight line segments that join the so-called Voronoi vertices 
(represented by asterisks). Each segment is closest to two exactly points and runs through part of 
the line that divides the plane at halfway between the two points. Some points are surrounded by 



®The algoritlmis as presented here are fairly inefficient — they run in O(n^) time and space. In practice, 0(n log n) 
times were achieved with appropriate data structures. 
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Algorithm 1 Prediction of midvein from MPEDs. Mode{i,j) denotes the action of PINl from cell 
i to cell j and can be: for no action, 1 if against gradient, and 2 if with gradient. 



1: 
2: 
3: 
4: 
5: 
6: 
7: 
8: 
9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26 



INPUT: Domain with MPED, n 
OUTPUT: Domain with MPED and midvein PED 
niidvcinincrcascd := true 
while midveinincreased == true do 
midveinincreased := false 
consistent := false 
while consistent == false do 
consistent := true 
Compute c as in Section 5 
for all Neighbors i and j do 

if c{i) > c{j) and Mode{i, j) ^ and Mode{i, j) ^ 1 then 

Mode{i,j) = 1, consistent :=false 
end if 

if c{i) < c{j) and Mode{i,j) ^ and Mode{i,j) ^ 2 then 

Mode{i,j) = 2, consistent :=false 
end if 
end for 
end while 

Find i,j neighbors and Ac = c{i) — c{j) where |Ac| is largest 
if Ac > Ti then 

Mode{i,j) = 1, midveinincreased=true 
end if 

if Ac < Ti then 

Mode{j,i) = 1, midveinincreased=true 
end if 
end while 
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Figure 21: How to use Voronoi diagrams to draw. (A) Points in 2-D. (B) Voronoi diagram of points 
in A. Voronoi vertices (asterisks) are joined by Voronoi edges or segments (lines). (C) Two groups 
of points: open circles and closed circles. Boundary between two groups (solid line) represents shape 
in D. (D) Shape from C. (E) Three groups of points. (F) Shape now consists of two regions that 
share an interface. The size of each region can be computed as the sum of sizes for each Voronoi 
cell, and the size of the interface is the sum of the sizes for each segment that comprises it. The 
Voronoi diagram and the convex hull of each Voronoi cell (volume and boundary) can be computed 
using the publicly available software package qhull or any of the algorithms in O'Rourke [40]. 
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these segments completely and thus define a closed region. This region, known as a Voronoi cell, is 
always convex and it is possible to determine its size (area). 

Now suppose that the points are separated into two groups: the open circles and the closed 
circles, as in Fig. 21 C. Then only some of the segments (solid lines) will separate the two groups, 
while the others will not (dashed lines). In this example, the points represented by a closed circle 
form a group that defines a closed region of space that consists of the union of the Voronoi cells 
of each of the points. This larger region (Fig. 21D) is no longer convex, but its size (area) can be 
computed by adding the sizes of the Voronoi cells. In addition, the size of the boundary of the 
region can be computed by adding the lengths of the segments that comprise it. Note that this 
procedure can be used to define any region of space, i.e. a shape, by using more points. 



Algorithm 2 ComputeSetup. The calculation of the Voronoi diagram as well as of the volume 
of the convex hull of a set of points is documented in O'Rourke [40]. The freely available software 
package qhull contains efficient implementations for both. 

1: INPUT: Groups of points Gi, . . . , G^, of which Gi is the bounding group; diffusion coefficient D 

2: OUTPUT: Geometry of fc — 1 cells CSegs, size S{i), production function p{i), interface size 
diffusion matrix M(i,j) 

3: { vsegs,vcells } ^ VORONOi(U*LiGj) 

4: for i = 2 to A; do 

5: S{i) ^ 

6: for all P G Gi do 

7: S{i) ■*— S(z)+CONVEXHULLVOLUME(t;ceWs(P)) 

8: end for 

9: CSegs(i - 1) ^ COMPUTEBouNDARYSEGS(Gi, tisegs,wce^ls) {See Algorithm 3} 
10: end for 

11: for z = 1 to fc — 1 do 

12: for j = 1 to fc — 1 and j ^ i do 

13: Hi J) ^ CoMPUTElNTERSECTiONLENGTH(Gi, Gj,i)ce//s) { See Algorithm 4} 

14: M{i,j) ^ D*I{i,j) 

15: end for 
16: end for 

17: for z = 1 to fc — 1 do 

18: M{l,l)^j: ^^-M{l,j) 

19: p{l) ^ 1/S{i) 

20: end for 



Algorithm 3 ComputeBoundarySegs(G, i^seps, vceZ/s). 

1: INPUT: A group of points G, Voronoi segments vsegs, Voronoi cells vcells 

2: OUTPUT: Segments segs that form the boundary of G with other groups. 

3: segs <— 

4: for seg e vsegs do 

5: if P G G,Q ^ G s.t. seg € vcell{P) and seg e vcell{Q) then 
6: segs <— segs U {seg} 

7: end if 

8: end for 
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Algorithm 4 COMPUTElNTERSECTIONLENGTH(Gi, Gj, wce^^s). 

1: INPUT: A group of points G, Voronoi segments vsegs, Voronoi cells vcells 
2: OUTPUT: Length of intersection len. 
3: len <— 

4: for seg G vsegs do 

5; if P G Gi,Q e Gj s.t. seg € vcell{P) and seg £ vcell{Q) then 
6: len ^ /en+LENGTH(ses') 

7: end if 
8; end for 

Now suppose that there are three groups of points as in Fig. 21E: open circles, closed circles, 
and closed diamonds. This configuration defines two closed regions of space (Fig. 21F) which share 
some of the Voronoi segments (dashed line). As before, the size of each region can be determined 
and the size of the interface between the two region can be computed as well. Thus, this method 
can not only define any shape but it can also provide both the shape size (area) and the size of 
the interface between any two neighboring shapes. These were the requirements for a collection of 
plant cells. 

Fig. 22 shows how this procedure is used to compute the geometry and topology of a traced 
root. The original image is of the root cell structure in a sUce of the Arabidopsis root hand-traced 
by an expert. In Fig. 22_B, we determine those pixels that belong to the interior of cells and color 
them white. Pixels with other colors separate different cells but are otherwise not used in the 
computation. We draw a bounding white curve around the root the pixels of which correspond to 
the points denoted by open circles in Fig. 21. We extract the coordinates of each white pixel and 
group the resulting points by computing the contiguous regions of white pixels. Thus, each cell 
interior is represented by a different (but unique) group of points and the bounding white curve is 
the only additional group. Then we follow the procedure above and obtain the shapes of each cell 
as well as their sizes and interfaces. The bounding curve is necessary in order to make sure that all 
Voronoi cells inside plant cells are finite regions of space. 

This same procedure works in higher dimensions as well. For example, the Voronoi diagram of 
points in three dimensions can be computed efficiently and three dimensional shapes (representing 
plant cell) can be obtained exactly similarly. The notion of cell size then becomes volume, and 
interface size becomes area. Both of these can be computed using the values from Voronoi cells 
as in the two-dimensional case discussed above. Thus, this approach produces the required graph 
structure for three-dimensional cell measurements as well. 
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Figure 22: Extracting the geometry and topology of a traced the root slice. (A) Image of a hand- 
traced root slice by an expert. (B) Processed image. White pixels in the interior of cells, other 
colors for display purposes and to separated different cells. (C) Result of procedure overlaid on 
top of original image. (D) Result of procedure. All cells are numbered but a dot "." is displayed 
for small cells. 
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